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The first stars fundamentally transformed the early universe by emitting the first 
light and by producing the first heavy elements. These effects were predetermined 
by the mass distribution of the first stars, which is thought to have been fixed by 
a complex interplay of gas accretion and protostellar radiation. We performed 
radiation-hydrodynamics simulations that followed the growth of a primordial 
protostar through to the early stages as a star with t her mo- nuclear burning. The 
circumstellar accretion disk was evaporated by ultraviolet radiation from the star 
when its mass was 43 times that of the Sun. Such massive primordial stars, in 
contrast to the often postulated extremely massive stars, may help explain the 
fact that there are no signatures of the pair-instability supernovae in abundance 
patterns of metal-poor stars in our galaxy. 

Theoretical studies and detailed computer simulations show that the cradles of the first stars were 
dense concentrations of primordial gas, with masses of ~ 1000 that of the sun. Such gas clouds formed 
through radiative cooling, with hydrogen molecules at the center of a dark matter halo of 10 6 solar mass 
(Mq), when the age of the universe was a few hundred million years old (1). 

According to our current understanding of star formation, a gas cloud's dense core gravitationally 
contracts in a non-homologous run-away fashion, in which the densest parts become denser faster than 
does the rest of the cloud {2-5). In a primordial gas cloud, one or a few embryo protostars are formed near 
the center (2, 6). The initial mass of these embryo protostars is only ~ 0.01 M©; the bulk of the dense 
core material remains in the surrounding envelope and is subsequently drawn toward the protostar (or 
protostars) through gravity. With the typical angular momentum of dense cores, the centrifugal barrier 
allows only a small amount of infalling gas to accrete directly onto the star. Instead, a circumstellar disk 
is formed and the gas is accreted onto the central star through the disk ( 7) . The final mass of these first 
stars is fixed, when the mass accretion terminates. However, when and how this termination occurs are 
largely unknown. 

Because the luminosity increases rapidly with protostellar mass, radiative feedback is expected to 
regulate the mass accretion and ultimately shut off the accretion flow, setting the final mass of the first 
stars. A primordial star more massive than a few tens times the mass of the Sun radiates a copious amount 
of hydrogen ionizing photons (> 13.6 eV). As an accreting star grows in mass, the ionized region in its 
vicinity grows, and eventually the circumstellar disk is directly exposed to the stellar ionizing radiation. 
The gas on the disk surface is photoionized and heated and evaporates away from the star-disk system. A 
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semi-analytical model of this process predicts that this photoevaporation quenches the accretion flow of 
the disk material and puts an end to the stellar growth (8). However, the interplay between the accretion 
flow and the stellar radiation is highly dynamical and complex. 

To identify the exact mechanism that halts the growth of a first star and to determine its final mass, 
we applied a method used for studying the present-day massive star formation (9, 10) to the case of the 
formation of the first stars in a proper cosmological context. We followed the radiation hydrodynamic 
evolution in the vicinity of an accreting protostar, incorporating thermal and chemical processes in 
the primordial gas in a direct manner. We also followed the evolution of the central protostar self- 
consistently by solving the detailed structure of the stellar interior with zero metallicity as well as the 
accretion flow near the stellar surface [supporting online material (SOM) text] (11, 12). We configured 
the initial conditions by using the results of a three-dimensional (3D) cosmological simulation, which 
followed the entire history from primeval density fluctuations to the birth of a small seed protostar at the 
cosmological redshift 14 (6). Specifically, when the maximum particle number density reached 10 6 cm~ 3 
in the cosmological simulation, we considered a gravitationally bound sphere of radius 0.3 pc around the 
density peak, which enclosed a total gas mass of ~ 300 Mq . We reduced the 3D data to an axisymmetric 
structure by averaging over azimuthal angles. 

The system was evolved until the central particle number density reached 10 12 cm~ 3 . We then 
introduced a sink cell of size ~ 10 astronomical units (AU) and followed the subsequent evolution of 
the central protostar with a stellar evolution code. The accretion rate onto the protostar was calculated 
directly from the mass influx through the sink-cell boundary, whereas the luminosity from the protostar, 
which controls the radiative feedback, was calculated consistently from the protostellar model by using 
the derived accretion rate. 

At its birth, a very small protostar of ~ 0.01 Mq was surrounded by a molecular gas envelope of 
~ 1 Mq, which was quickly accreted onto the protostar. Atomic gas further out initially had too much 
angular momentum to be accreted directly, and a circumstellar disk formed. The infalling atomic gas first 
hit the disk plane roughly vertically at supersonic velocities. A shock front formed; behind the shock, the 
gas cooled and settled onto the disk, and its hydrogen was converted to the molecular form via rapid gas 
phase three-body reactions. The molecular disk extended out to ~ 400 AU from the protostar, when the 
stellar mass was 10 Mq. Accretion onto the protostar proceeded through this molecular disk as angular 
momentum was transported outward. The accretion rate onto the protostar was ~ 1.6 x 10~ 3 M Q year" 1 
at that moment. 

The evolution of the central protostar is determined by competition between mass growth by accretion 
and radiative energy loss from the stellar interior. The time scale for the former is the accretion timescale 
t acc = M*/M, where M* is mass of the protostar and M is mass accretion rate, whereas that for the 
latter is the Kelvin- Hclmholtz (KH) timescale £kh = GAf 2 /i?*L*, where L* is the luminosity from the 
stellar interior, i?„ is radius of the star, and G is gravitational constant. The total luminosity of the 
protostar L to t is the sum of the stellar luminosity and accretion luminosity L acc = GM*M / Ft* (13). 

In the early accretion phase, the stellar radius remained almost constant at ~ 50 solar radius (R®) 
(Fig. [T]A). The stellar luminosity was substantially lower than the accretion luminosity (Fig. QjB) and 
the KH timescale was much longer than the accretion time (Fig. [U-C). Consequently, entropy carried by 
the accreted gas accumulated at the stellar surface nearly without loss. During this quasi-adiabatic stage 
(M* < 7 Mq), the luminosity L* increased with stellar mass. When the star grew to 8 Mq, the KH 
timescale finally fell below the accretion timescale (Fig. [TJC). After this, the protostar began its so-called 
KH contraction, in which it gradually contracted as it radiated its energy away (Fig. Q3A). The stellar 
luminosity was the main component of the total luminosity after this evolutionary stage (Fig. [T]-B). The 
stellar luminosity increased, and stellar radius i?* decreased, as the stellar mass increased. As a result, 
the stellar effective temperature T e g oc (L^/i? 2 ) 1 / 4 and the ultraviolet (UV) flux rapidly rose (Fig. [1}D). 
Thus, ionization and heating by UV photons became important already in the KH contraction stage. 

When the stellar mass was 20 Mq , an ionized region rapidly expanded in a bipolar shape perpendicular 
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to the disk, where gas was cleared away (Fig. H|A). At this moment, the disk extended out to ~ 600 AU. 
The disk was self-shielded against the stellar H2-dissociating (11.2 eV < hv < 13.6 eV) as well as the 
ionizing radiation. The ionized atomic hydrogen (HII) region continued to grow and finally broke out 
of the accreting envelope. At A/* ~ 25 Af©, the size of the bipolar HII region exceeded 0.1 pc (Fig. 
I2J-B). Because of the high pressure of the heated ionized gas, the opening angle of the ionized region 
also increased as the star grew (Fig. [2]-C). Shocks propagated into the envelope preceding the expansion 
of the ionized region. The shocked gas was accelerated outward at a velocity of several kilometers per 
second. The shock even reached regions shielded against direct stellar UV irradiation. The outflowing 
gas stopped the infall of material from the envelope onto the disk (Fig. S6). Without the replenishment 
of disk material from the envelope the accretion rate onto the protostar decreased (Fig. [3]) . In addition, 
the absence of accreting material onto the circumstellar disk means that the disk was exposed to the 
intense ionizing radiation from the star. The resulting photoevaporation of disk gas also reduced the 
accretion rate onto the protostar. The photoevaporated gas escaped toward the polar direction within 
the ionized region. The typical velocity of the flow was several tens of kilometers per second, comparable 
with the sound speed of the ionized gas, which was high enough for the evaporating flow to escape from 
the gravitational potential well of the dark matter halo. 

When central nuclear hydrogen burning first commenced at a stellar mass of 35 Af©, it was via the 
proton-proton (pp) -chain normally associated with low-mass stars. The primordial material does not 
have the nuclear catalysts necessary for carbon-nitrogen-oxygen (CNO) -cycle hydrogen burning. Because 
the pp-chain alone cannot produce nuclear energy at the rate necessary to cover the radiative energy loss 
from the stellar surface, the star continues to contract until central temperatures and densities attain 
values that enable the 3-a process of helium burning (11). The product of helium burning is carbon, 
and once the relative mass abundance of carbon reaches ~ 10 -12 , CNO-cycle hydrogen burning takes 
over as the principal source of nuclear energy production, albeit at much higher central densities and 
temperatures than those of stars with solar abundances. These first-generation ZAMS (zero-age main 
sequence) stars are thus more compact and hotter than their present-day counterparts of equal mass (14). 
The subsequent evolution of the accreting star followed along the ZAMS mass-radius relationship (Figure 
UK). By the time the star attained 40 Af©, the entire region above and below the disk (Fig. [2jD) was 
ionized. Mass accretion was terminated when the stellar mass was 43 Af© (Fig. [3]). 

The entire evolution described above took about 0.1 million years from the birth of the embryo 
protostar to the termination of the accretion. The star is expected to live another few million years 
before exhausting all available nuclear fuel and exploding as a core-collapse supernova (15). 

Our calculations show that the first stars regulated their growth by their own radiation. They were 
not extremely massive, but rather similar in mass to the O-type stars in our Galaxy. This resolves a 
long-standing enigma regarding the elemental abundance patterns of the Galactic oldest metal-poor stars, 
which contain nucleosynthetic signatures from the earliest generation of stars. If a substantial number of 
first stars had masses in excess of 100 Af©, they would end their lives through pair- instability supernovae 
(16, 17), expelling heavy elements that would imprint a characteristic nucleosynthetic signature to the 
elemental abundances in metal-poor stars. However, no such signatures have been detected in the metal- 
poor stars in the Galactic halo (18, 19). For example, the abundances of elements with odd atomic 
numbers are generally reduced in remnants of primordial supernovae. The odd-even contrast pattern 
expected in pair-instability supernovae is much stronger than the observed patterns in Galactic metal- 
poor stars (17). Moreover, pair-instability supernovae predict a small abundance ratio [Zn/Fe], but 
observed values are much larger (16). Detailed spectroscopic studies of extremely metal-deficient stars 
indicate that the metal-poor stars were born in an interstellar medium that had been metal-enriched by 
supernovae of ordinary massive stars (20). 

Second-generation stars, which formed from the primordial gas affected by radiative or mechanical 
feedback from the first stars, could have dominated the metallicity of the young interstellar medium, 
which then spawned the observed Galactic halo stars. These second-generation stars could have been 
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more numerous but less massive than the first stars because of a different gas thermal evolution, with 
additional radiative cooling via H 2 and HD molecules (21, 22). However, this mode of star formation 
is suppressed even with weak H 2 photodissociating background radiation (22). If so suppressed, the 
formation process of the later-generation primordial stars would be similar to that of the very first stars, 
and most primordial stars could have experienced the evolution presented in this article regardless of 
their generation. One might argue that today's observed metal-deficient stars formed after an episode 
of star formation with non-zero metallicity. For the star formation process to differ substantially from 
that of the very first stars, one would require metallicities in excess of [Fe/H] > —5 ( 23, 24). Caffau 
and co-workers (25) report on observations of a metal-deficient star with [Fe/H] ~ —5 but without 
the corresponding enhancement of carbon, nitrogen, and oxygen found in metal-deficient stars. The 
abundance pattern of this star agrees with expectations from core-collapse supernovae, implying that it 
formed from gas enhanced by material ejected from primordial stars with masses less than 100 M©. 

We have performed radiation-hydrodynamic simulations only for a single star-forming region embed- 
ded in a cosmological simulation. Our selected dark halo was typical in mass, spin, and formation epoch, 
when compared with those in other studies (5). The evolution presented here is somewhat similar to 
that predicted by the semi-analytic model, in which the expansion of the ionized region begins when the 
stellar mass is ~ 25 M Q and the final mass is ~ 57 M© (8), the lowest final stellar mass predicted by 
the semi-analytic treatment. If input parameters of the semi-analytic models are chosen to fit our initial 
gas cloud, however, the final mass should be higher, ~ 90 M© (SOM text). Our calculations follow the 
dynamical response of the infalling gas onto the circumstellar disk. The expansion of the ionized region 
around the protostar generates a powerful outflow even behind the surrounding disk. This effect reduces 
the accretion rate significantly. Thus, our radiation-hydrodynamic calculations predict systematically 
lower final masses than those of the semi-analytic models. 

Although the results described above provide a complete picture of how a primordial protostar regu- 
lates and terminates its growth, there are a few key quantities that determine the strength of the feedback 
effect, as suggested by the semi-analytic model (8). With smaller initial rotation of the natal dense core, 
the density in the envelope near the polar directions would be higher. Hydrogen recombination occurs 
rapidly in the dense gas, which prevents the breakout of the ionized region. Gas accretion can last for a 
longer time in this case, and would form more massive stars (SOM text). Nevertheless, even considering 
variations among dark halos bearing the first stars, a substantial fraction of the first stars should be 
less massive than 100 M© and end their lives as ordinary core-collapse supernovae. Gas accretion might 
not be completely halted in a few exceptional cases, thus leading to the formation of a small number of 
extremely massive stars that are > 100 M© in the early universe (26). Black holes left after such very 
massive stars' deaths might have grown up to be the supermassive black holes lying in galaxies. 

Recent 3D cosmological simulations showed that a primordial gas cloud breaks up into several embryo 
protostars in an early phase (27, 28). Each of these protostars would continue to grow through mass 
accretion, but it remained uncertain how and when the growth is halted. Our radiation-hydrodynamics 
calculations explicitly show that the parent gas cloud is evaporated by intense radiation from the cen- 
tral star when its mass is several tens of solar-masses. The circumstellar disks in our simulations are 
marginally stable against gravitational fragmentation [Fig. S5] . We expect that these disks - in a 3D sim- 
ulation - would evolve analogously to our numerical simulations with assumed axial symmetry and have 
a similar time-averaged structure. It is conceivable that a few protostars are ejected dynamically from 
the parent cloud, to remain as low-mass stars. Observationally, however, there have been no low-mass 
zero-metallicity stars discovered in the Galaxy. This fact suggests the limited formation efficiency of such 
low-mass primordial stars (19). Low-mass stars (< 1 M©) and extremely high mass stars (> 100 M©), if 
any, are thus a minor population among the first stars. 

Our self-consistent calculations show that the characteristic mass of the first stars is several tens of 
solar masses. Although this is less than that of the conventionally proclaimed extremely massive stars 
(> 100 M©), it is still much larger than the characteristic mass of stars in our galaxy (< 1 M©) (29). 



4 



References and Notes 

1. V. Bromm, N. Yoshida, L. Hcrnquist, C. F. McKee, Nature 459, 49 (2009). 

2. K. Omukai, R. Nishi, Astrophys. J. 508, 141 (1998). 

3. T. Abel, G. L. Bryan, M. L. Norman, Science 295, 93 (2002). 

4. N. Yoshida, K. Omukai, L. Hernquist, T. Abel, Astrophys. J. 652, 6 (2006). 

5. B. W. O'Shca, M. L. Norman, Astrophys. J. 654, 66 (2007). 

6. N. Yoshida, K. Omukai, L. Hernquist, Science 321, 669 (2008). 

7. H. W. Yorkc, P. Bodenheimer, Astrophys. J. 525, 330 (1999). 

8. C. F. McKee, J. C. Tan, Astrophys. J. 681, 771 (2008). 

9. H. W. Yorke, A. Welz, Astron. & Astrophys. 315, 555 (1996). 

10. H. W. Yorkc, C. Sonnhalter, Astrophys. J. 569, 846 (2002). 

11. K. Omukai, F. Palla, Astrophys. J. 589, 677 (2003). 

12. T. Hosokawa, K. Omukai, Astrophys. J. 691, 823 (2009). 

13. This definition of accretion luminosity includes the mechanical luminosity of an accretion-driven wind 
and the ratio tua/tucc is equal to the ratio L acc /i» . 

14. D. Ezcr, A. G. W. Cameron, Astrophys. & Space Science 14, 399 (1971). 

15. D. Schaercr, Astron. & Astrophys. 382, 28 (2002). 

16. H. Umeda, K. Nomoto, Astrophys. J. 565, 385 (2002). 

17. A. Hcgcr, S. E. Woosley, Astrophys. J. 567, 532 (2002). 

18. J. Tumlinson, A. Venkatesan, J. M. Shull, Astrophys. J. 612, 602 (2004). 

19. A. Frebcl, J. L. Johnson, V. Bromm, Mon. Not. R. Astron. Soc. 392, L50 (2009). 

20. N. Iwamoto, H. Umeda, N. Tominaga, K. Nomoto, K. Macda, Science 309, 451 (2005). 

21. B. W. O'Shea, T. Abel, D. Whalen, M. L. Norman, Astrophys. J. Lett. 628, L5 (2005). 

22. N. Yoshida, K. Omukai, L. Hernquist, Astrophys. J. Lett. 667, L117 (2007). 

23. [Fe/H] = log(Fe/H) sta r - log(Fc/H) sun , where (Fe/H) is mass ratio of Fc to H . 

24. T. Hosokawa, K. Omukai, Astrophys. J. 703, 1810 (2009). 

25. E. Caffau, et al, Nature 477, 67 (2011). 

26. T. Ohkubo, K. Nomoto, H. Umeda, N. Yoshida, S. Tsuruta, Astrophys. J. 706, 1184 (2009). 

27. M. J. Turk, T. Abel, B. O'Shea, Science 325, 601 (2009). 

28. P. C. Clark, et al, Science 331, 1040 (2011). 

29. G. Chabrier, Pub. Astron. Soc. Pac. 115, 763 (2003). 

5 



30. S. Richling, H. W. Yorke, Astrophys. J. 539, 258 (2000). 

31. H. W. Yorke, M. Kaisig, Computer Physics Communications 89, 29 (1995). 

32. D. E. Osterbrock, G. J. Ferland, Astrophysics of gaseous nebulae and active galactic nuclei (2006). 

33. D. Hollenbach, C. F. McKee, Astrophys. J. Suppl. Ser. 41, 555 (1979). 

34. P. R. Shapiro, H. Kang, Astrophys. J. 318, 32 (1987). 

35. P. Anninos, Y. Zhang, T. Abel, M. L. Norman, New Astronomy 2, 209 (1997). 

36. C. D. Levermore, G. C. Pomraning, Astrophys. J. 248, 321 (1981). 

37. N. J. Turner, J. M. Stone, Astrophys. J. Suppl. Ser. 135, 95 (2001). 

38. M. Mayer, W. J. Duschl, Mon. Not. R. Astron. Soc. 358, 614 (2005). 

39. B. T. Draine, F. Bcrtoldi, Astrophys. J. 468, 269 (1996). 

40. J. C. Tan, C. F. McKee, Astrophys. J. 603, 383 (2004). 

41. T. Abel, P. Anninos, Y. Zhang, M. L. Norman, New Astronomy 2, 181 (1997). 

42. D. Galli, F. Palla, Astron. & Astrophys. 335, 403 (1998). 

43. F. Palla, E. E. Salpeter, S. W. Stahler, Astrophys. J. 271, 632 (1983). 

44. A. G. G. M. Tielens, D. Hollenbach, Astrophys. J. 291, 722 (1985). 

45. A. Stacy, T. H. Grcif, V. Bromm, Mon. Not. R. Astron. Soc. 403, 45 (2010). 

46. N. I. Shakura, R. A. Sunyaev, Astron. & Astrophys. 24, 337 (1973). 

47. M. R. Krumholz, R. I. Klein, C. F. McKee, S. S. R. Offner, A. J. Cunningham, Science 323, 754 
(2009). 

48. T. Hosokawa, H. W. Yorke, K. Omukai, Astrophys. J. 721, 478 (2010). 

49. S. W. Stahler, F. H. Shu, R. E. Taam, Astrophys. J. 241, 637 (1980). 

50. M. N. Machida, K. Omukai, T. Matsumoto, S.-i. Inutsuka, Astrophys. J. Lett. 647, LI (2006). 

51. J. C. Tan, E. G. Blackman, Astrophys. J. 603, 401 (2004). 

52. J. Silk, M. Langer, Mon. Not. R. Astron. Soc. 371, 444 (2006). 

53. J. C. Tan, B. D. Smith, B. W. O'Shea, American Institute of Physics Conference Series, D. J. Whalen, 
V. Bromm, & N. Yoshida, ed. (2010), vol. 1294 of American Institute of Physics Conference Series, 
pp. 34-40. 

54. M. N. Machida, S.-i. Inutsuka, T. Matsumoto, Astrophys. J. 724, 1006 (2010). 

55. K. M. Kratter, C. D. Matzner, M. R. Krumholz, R. I. Klein, Astrophys. J. 708, 1585 (2010). 

56. D. Hollenbach, D. Johnstone, S. Lizano, F. Shu, Astrophys. J. 428, 654 (1994). 

57. F. D. Kahn, Astron. & Astrophys. 37, 149 (1974). 

58. H. W. Yorke, E. Kruegel, Astron. & Astrophys. 54, 183 (1977). 

6 



59. M. G. Wolfire, J. P. Cassinelli, Astrophys. J. 319, 850 (1987). 

60. R. Kuipcr, H. Klahr, H. Bcuthcr, T. Henning, Astrophys. J. 722, 1556 (2010). 

61. R. Kuiper, H. Klahr, H. Bcuther, T. Henning, Astrophys. J. 732, 20 (2011). 

62. T. Peters, et al, Astrophys. J. 711, 1017 (2010). 

63. J. E. Dale, I. A. Bonnell, C. J. Clarke, M. R. Bate, Mon. Not. R. Astron. Soc. 358, 291 (2005). 

64. T. Peters, R. S. Klcssen, M.-M. Mac Low, R. Banerjee, Astrophys. J. 725, 134 (2010). 

Acknowledgments: We thank T. Nakamura, K. Nomoto, S. Inutsuka, and N. Turner for stimulat- 
ing discussions on this topic. Comments by an anonymous referee helped improve the manuscript. T.H. 
appreciates the support by Fellowship of the Japan Society for the Promotion of Science for Research 
Abroad. The present work is supported in part by the grants-in-aid by the Ministry of Education, Science 
and Culture of Japan (19047004, 2168407, 21244021:KO, 20674003:NY). Portions of this research were 
conducted at the Jet Propulsion Laboratory, California Institute of Technology, which is supported by 
NASA. Data analysis was (in part) carried out on the general-purpose PC farm at Center for Computa- 
tional Astrophysics (CfCA) of National Astronomical Observatory of Japan. 



Supporting Online Material 

www . sciencemag . org 
SOM text 
Figures S1-S7 
Tables SI and S2 
References (30-64) 

25 April 2011; accepted 28 October 2011 
Published online 10 November 2011 



7 



(a) Radius 




(b) Bolometric luminosity 
■ ' 



(c) Timescales 




(d) UV luminosities 




ionizing 



10 30 
Stellar mass: log M t ( M s ) 



Fig. 1: Evolution of the stellar radius, bolometric luminosity, evolutionary timescales, and ionizing 
(hv > 13.6 eV) and dissociating (11 eV < hv < 13.6 eV) photon number luminosity. In panel (b), 
the total luminosity L tot (black) is the sum of the stellar luminosity (red) and accretion luminosity 
Luce (blue). The KH timescale iKH (red) and accretion timescale t acc (blue) are depicted in panel (c). 
The yellow and blue backgrounds denote the adiabatic accretion phase and KH contraction phase in the 
protostellar evolution. 
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Fig. 2: UV radiative feedback from the primordial protostar. The spatial distributions of gas temperature 
(left), number density (right), and velocity (right, arrows) are presented in each panel for the central 
regions of the computational domain. The four panels show snapshots at times, when the stellar mass 
is M* = 20 M Q (a), 25 M (b), 35 M (c), and 42 M (d). The elapsed time since the birth of the 
primordial protostar is labeled in each panel. 
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10 20 30 40 50 60 

Stellar Mass: M*(M© ) 

Fig. 3: Evolution of the accretion rate onto the primordial protostar. The blue line indicates the evolution, 
which includes the effect of UV radiative feedback from the protostar. The red line indicates a numerical 
experiment with no UV feedback. The open and solid circles denote the characteristic epochs of the 
protostellar evolution, beginning of the KH contraction and the protostar's arrival to the ZAMS. Figure 
[2J A to D, shows the snapshots at the moments marked here with asterisks. 
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Supporting Online Material 



1 Numerical Methods 

We have developed a hybrid code to simultaneously calculate the evolution of an accreting first generation 
(proto-)star and the flow of the surrounding primordial material. We solve the accretion flow with a 
radiation hydrodynamic code, whereby the central (proto-)star is represented as a sink cell. This sink cell 
grows in mass by inflow through its boundary and it influences the surrounding material flow via gravity 
and by emitting radiation. We follow the evolution of the central star from the earliest protostellar stage 
up to the early phases of nuclear burning on the zero age main sequence by solving the conventional 
stellar structure equations, taking into account mass accretion. The material functions for the stellar 
evolution calculation - the opacity, the equation of state and other thermodynamic relationships, as well 
as the nuclear reaction network - assume zero metallicity. The calculated evolution provides the emitted 
radiation from the sink cell as a function of time. We describe below our methods of calculation for the 
accretion flow and the protostar, respectively. 



1.1 Radiation Hydrodynamics of the Accretion Flow 

To study the evolution of the flow of primordial gas we employ a grid-based axisymmetric radiation 
hydrodynamic code with self-gravity, previously used for studying present-day star formation (7,10) and 
the evolution of photoionized gas flow from protostellar disks (9, 30). This code makes use of a nested-grid 
technique to cover a wide spatial range (31). We have added a chemical network besides changes based 
on the nature of the primordial material. 

The governing equations for gas dynamics in cylindrical coordinate (R, Z) are, 

g + v-(H = o, (i) 

d(ov) _ _ A 2 

-g-^ + V • (pv ® v) = -pV$ - Vp + n R + K (2) 

3c 

— + V • (ev) = -pV • v + Y - A, (3) 

p=(7-l)e, (4) 

where p and p are the gas density and pressure, $ the gravitational potential, v the 2-dimensional (2D) 
velocity vector (vr,Vz), A = pRv^ the angular momentum per unit volume, the radial unit vector, 
K the radiation force, e the gas internal energy density, Y and A the heating and cooling rates per unit 
volume, and 7 the adiabatic exponent. We calculate 7 from the chemical composition at each grid cell as 
in (2). The radiative and chemical processes included in the energy source term T — A are summarized 
in Table SI. 



1.1.1 Radiation Transfer 

For the problem at hand photons in different wavelength regimes play vastly different roles in the ther- 
mal and chemical processes. For example, whereas gas cools via molecular hydrogen line emission and 
infrared/optical continuum radiation, photons with energies hv > 13.6 eV (EUV) photoionize atomic 
hydrogen and photons with energies 11.2 eV < hv < 13.6 eV (FUV) photodissociate molecular hydrogen. 
To treat their transfer properly without time-consuming calculations, we adopt different methods for each 
of the radiation components. 
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Table SI. Included thermal processes 





Processes 


References 




... 

Photoionization 


(32) 


Heating 


Photodissociation 


(33) 




H2 formation 






H2 collisional excitation 


§U.1.1| 




H - free-bound emission 


§11-1-11 




H2 collisional dissociation 


(33), (34) 


Cooling 


H collisional ionization 


(32) 




H collisional excitation 


(35) 




Compton scattering 


(35) 




Hell collisional excitation 


(35) 



Molecular Hydrogen Line Cooling The primary cooling process in low-temperature (<~ 8000K) 
primordial gas is line emission via rotational and vibrational transitions of hydrogen molecules. We adopt 
the fitting formula by (33) for the optically-thin limit for n <~ 10 9 cm~ 3 . Dense gas with n >~ 10 9 cm' 3 
is opaque to H2 line emission, and the cooling rate is reduced by photon trapping. In this case the cooling 
rate is calculated by adapting the method of (4) to the case of axial symmetry. We sum up cooling rates 
for all possible transitions among the rotational levels from J = to 20 and vibrational levels v = 0,1, 2. 
The cooling rate by a transition for the optically thick case is obtained by multiplying the value for 
optically thin cooling by the escape probability. The escape probability [3 esc is evaluated as 



P(tr) + P(tz) 



(5) 



where /3(t) is 



1 — exp(r) 



(6) 



and tr and tz are the optical depths of this transition along R and Z directions, respectively. We 
calculate the optical depths using the local velocity gradients (so-called Sobolev approximation) 

(7) 



'\dv q /dq\' 



where q = R or Z , a is the absorption coefficient and c s is the sound speed. 

For calculating the level-population of H2 molecules, we omit the excitation by absorbing stellar 
photons, which could potentially reduce the line cooling rate. In general, H2 molecules could be excited 
to (i) the Lyman- Werner bands by absorbing FUV photons or (ii) the higher rotational and vibrational 
levels by absorbing infrared photons. In our case, the process (i) would be minor because the disk is 
shielded against the stellar FUV radiation (also see Sec. The process (ii) is negligible in most parts 
of the disk even with the radiation from the protostar. 



Continuum Cooling During the collapse phase before the formation of a protostar, continuum cooling 
via H2 collision-induced emission (CIE) is important in the dense (n > 10 13 cm -3 ) molecular gas (2). 
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However, we only find such dense gas inside the sink cell; it does not otherwise appear in our hydro- 
dynamical calculation. Instead, cooling via H - free-bound emission becomes important in the nearly 
vertical flows onto the circumstellar disk (also see Sec. |3] below). 

The continuum cooling rate A c in the gray approximation is written as 

A c = cpn P (^B(T g ) - E c ), (8) 

where c is the speed of light, and up is the Planck mean opacity per unit mass. The continuum radiation 
energy density E c satisfies the equation 

dE 

= -V ■F c + A c +j, 1 (9) 

where F c is the energy flux, and j* is the stellar source term. Using the total luminosity of the protostar 
Ltot, the stellar source term is given by 

J * = 2nAR?AZ' (10) 

where AR and AZ are the size of the sink cell, which is equal to the cell size at the finest grid-level, 
12 AU. The source term is zero except for the sink cell. We solve equations (J9j) with the flux-limited 
diffusion (FLD) approximation (36) using the operator splitting technique with de/dt = — A c (37). The 
FLD method adopts a closure relation between F c and E c 

F c = - — VE C , (11) 

PUR 

where kr is the Rosseland mean opacity, and A is the flux-limiter defined as 

A = 6T5^ < 12 » 

E c pn R 

We use the opacities Kp and kr for primordial gas calculated by (38) in tabulated forms. As mentioned 
above, the most important continuum cooling is that via H~ free-bound emission. Neutral gas falling 
onto the disk is heated up to T g ~ 5 x 10 3 K by compressional heating, until thermal balance is achieved 
with this cooling process. In our calculation the accreting envelope is always optically thin to the non-UV 
continuum radiation and non-UV stellar radiation escapes without significant heating in the envelope. 



Stellar EUV and FUV Radiation EUV radiation from the star ionizes the material in its immediate 
vicinity. Within the HII region, there are two components of the EUV field: the direct stellar component 
and the diffuse component emitted via recombinations directly into the ground state (i.e. no "on-the- 
spot" approximation). We solve the transfer of these two EUV components separately following (30). 
We adopt a frequency-averaged approximation for each component, taking into account the difference in 
their mean energies, /iP* (stellar) and hvd (diffuse). 

To calculate the direct EUV field we cast a number of radial rays from the central star to the outer 
edge of the simulation box, along which the transfer equation for the direct EUV photon number flux IF* 
is solved: 

V r • ZF* = —n(l — x)ct*^"», (14) 

where x is the degree of ionization and cr» is the absorption cross section per particle. The cross section 
it* is a function of mean energy of the direct EUV photons, i.e., the effective temperature of the star 
T cff (32). 
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We employ the FLD approximation for the diffuse component, whose transfer equation is 



<>Ad =-V-F d + ai (T g )n 2 x 2 - n(l - x)a d chf d , (15) 



dt 



where N d is the diffuse EUV photon number density, a± is the recombination coefficient to the ground 
state, and <r d is the absorption cross section. Under the approximation that the mean energy of the 
diffuse component depends only on the local gas temperature T g , the cross section a d becomes a function 
of T g (32). In the FLD approximation, 

T d = - CX VAf d , (16) 
n(l - x)(j d 

in analogy to equation (fTTj) . The photoionization heating and recombination cooling rates are calculated 
for each component of EUV radiation separately. 

As for the photodissociating FUV radiation, we only consider the stellar direct component. We 
calculate the photodissociation rate by multiplying the value for the optically-thin case by the self- 
shielding factor (39) 

r i (^h 2 < Nx) 

fsh(Nn 2 ) = { (NnjN,)- 3 / 4 (JVi < N H , < N 2 ) (17) 
{ (N H . 2 > N 2 ), 

where Nn 2 is the H2 column density from the protostar to the point under consideration, N\ = 10 14 cm -2 , 
and N2 = 10 22 cm -2 . Although the original functional form of f a h(N-a 2 ) by (39) is for the range 
N\ < Ah 2 < N2, we adopt the cut-off at Nn 2 > N2, where the photodissociation rate falls sharply 
according to (39). In our calculation, the column density largely exceeds N2 along the radial rays 
incident on the circumstellar disk. 

As shown in Sec.[3]below, the innermost part of the disk at R < several 10 AU is not spatially resolved 
in our calculations. We expect the innermost part of the disk to be optically thick and shield the material 
near the disk mid-plane behind it from the stellar radiation (40). We adopt the following treatment for 
modeling this effect. During the calculation, we evaluate the disk scale height H(R) = c s /Q at each 
radius R on the equator. We compute in,d, the i?-index of the cell at the outer edge of the unresolved 
part, for which H(R) < AZ. The unresolved part is assumed to be opaque to radial rays incident on the 
(iR,d + l)-th cell on the equator up to a height AZ, but transparent for the other rays. 



Radiation Force The total radiation force K is calculated as the sum of contributions by the contin- 
uum radiation, direct and diffuse EUV radiation: 

pn R n(l - x)o* n(l - x)g d 
K = F c H hv^F* H hv d F d . (18) 

c c c 

Here, we omit contributions from FUV photons. In reality, the gas accretion envelope is opaque against 
Lyman-a photons emitted from the stellar atmosphere and HII region. Radiation pressure via the Lyman- 
a scattering is exerted on the accretion envelope. With mass accretion through the circumstellar disk, in 
particular, photons are preferentially transferred toward the polar direction, where the gas density rapidly 
decreases as the mass accretion proceeds ("flashlight effect": c.f. (7)). The semi-analytic modeling by (8) 
shows that the Lyman-a pressure could influence the dynamics of infalling material near the rotation axis. 
As also discussed in (8), however, the radiation pressure via Lyman-a scattering would be significantly 
reduced once gas in polar directions is blown away and photons escape from the cavity. In this paper, 
we focus on the UV radiative feedback effects to derive upper limits of the stellar final masses. 
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Table S2. Included chemical reactions 



No. 



Reactions 



References 



Rl 

R2 

R3 

R4 

R5 

R6 

R7 

R8 

R9 

RIO 

Rll 

R12 

R13 

R14 

R15 



H + c^H++2c 
H+ + e -> H + 7 
H" + H -> H 2 + c 



(P) 
(32) 

(42) 
(42) 
(42) 
(34) 
(43) 
(43) 
(43) 
(42) 
(43) 
(32) 



H 2 + H+ -> H+ + H 



H 2 +e->2H|e 
H 2 + H — > 3 H 



3 H -> H 2 + H 



2 H + H 2 -> 2 H 2 
2 H 2 — s> 2 H + H 2 
H + c -> H~ + 7 
2H^H++e + H 
H + 7 — > H+ + e 
H 2 + 7 -> 2 H 



(44) 



and (fTfjl 



H-+c^H + 2e 
H" + H+ -> 2 H 



(^) 
(^) 



1.1.2 Chemical Reactions in the Primordial Gas 

We consider the five species of primordial gas H, H + , e, H 2 , and H~, based on the minimal model of (41) 
with some additional reactions. The included reactions are summarized in Table S2. Except H~, we 
solve kinetic equations with an implicit difference scheme. H - is assumed to be in chemical equilibrium. 
Photodissociation of H~ is omitted in the adopted chemical network. In our calculations the gas density is 
high enough that the collisional destruction processes (R4, 14, and 15) control its abundance. We curtail 
helium chemistry by assuming that helium is singly ionized in an HII region and atomic elsewhere. We 
do not include deuterium reactions because HD cooling is relevant only in low-temperature (T g < 200 K) 
and low-density (n < 10 8 cm -3 ) gas, which does not appear in our simulations. 

1.1.3 Angular Momentum Transport in the Accretion Disk 

In an exactly axisymmetric system, angular momentum must be conserved. In reality, however, a massive 
circumstellar disk can develop a non-axisymmetric spiral pattern, which exerts a torque on the matter 
in the disk and transfers angular momentum outward. Recent 3-dimensional (3D) numerical simulations 
demonstrated that this mechanism operates in fact in the circumstellar disks of the first stars (28,45). 
We mimic this effect by adopting the angular momentum transport via the so-called a- viscosity (46). 
The equation of angular momentum transport is thus given by 




(19) 



where Q is angular velocity, 77 = 2apc 2 s / (30), and a is a dimensionless free parameter. 
We assume that the a-parameter depends on the height from the equator, 




(20) 
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where ao is a constant. The effective values of the a-parameter in rapidly accreting circumstellar disks 
are a ~ 0.1 — 1 as estimated from 3D simulations of present-day massive star formation (4T) as well as the 
formation of the first stars (28). In the fiducial case explained in the main article, we adopted ao = 0.6. 

Figure [ST] shows the evolution of accretion rates onto the protostar for different values of ao- Although 
the evolution is qualitatively similar in all cases, with higher ao angular momentum is transferred more 
rapidly and therefore accretion rates onto the protostar become higher. We see that the final stellar mass 
increases with ao- With the low value of ao = 0.3, the final stellar mass is about 35 M Q . This is close 
to the low-mass limit of 30 M Q predicted by (3), who expected that this amount of gas would accrete 
onto the protostar in a few thousand years, which is much shorter than the timcscalc of the protostellar 
evolution. Even with ao as large as unity, on the other hand, the final mass reaches at most ~ 50 M@. 
Note, however, that this does not exclude the possibility of formation of higher-mass stars. The semi- 
analytic models predict that more massive stars would form with weaker initial rotation of the natal core 
(8) (also see Sec. [3]). Although rare, stars exceeding 100 M Q might still form in such circumstances. 

Figure IS 1 1 also shows the evolution in such a test case, whereby the initial angular momentum is 
artificially reduced to 30% of the original valued The final stellar mass is ~ 85 M Q in this case. The 
increase of the final mass is understood as follows. First, the gas density remains high in the accretion 
envelope in the polar directions for the same stellar mass. The stellar EUV photons are consumed more 
effectively by photoionization of neutral hydrogen generated via rapid recombination, which delays growth 
of an HII region. Second, the protostellar evolution differs from the fiducial case, because of the higher 
accretion rates. Figure [S2l shows that, with the weaker rotation, the Kelvin-Helmholtz contraction stage 
is shifted to higher stellar masses. Correspondingly, the stellar EUV luminosity increases rapidly at a 
stage when the protostar is more massive than for the fiducial case. Because of these effects, the stellar 
UV feedback becomes effective at higher stellar masses. We note that for all examined cases, the mass 
accretion ceases soon after the protostar's arrival to the zero age main sequence. 



1.2 Protostellar Evolution 

We follow the evolution of the central protostar by solving the four stellar structure equations taking 
account of mass accretion (11, 12,48): 



( dr\ _ 1 



\dMJ t Anpr 2 ' 



(21) 



(8P\ __GM 
[0L\ =e _ T f9_s 



where M is the Lagrangian mass coordinate, e is the energy production rate by nuclear fusion, s is the 
specific entropy, and L s is the radiative luminosity with adiabatic temperature gradient. The coefficient, 
C in equation is unity if L < L s (i.e., in radiative layers), and given by the mixing- length theory if 
L > L s (i.e., in convective layers). We also solve the structure of the accretion flow inside the sink cell 
under the assumption of steady state and spherical symmetry. The entire structure of both the protostar 
and accretion flow is consistently determined to satisfy the jump conditions for the accretion shock at the 
stellar surface (49). Mass accretion rates onto the protostar M* are given by mass inflow rates through 



lr The angular momentum is reduced at the beginning of the 2D calculation, e.g., at a point when the central density is 
10 6 cm -3 in the run-away collapse stage (see Sec.[2]|. 
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the surface of the sink cell in the radiation-hydrodynamics calculation. With a high accretion rate of 
M* > Mq yr _1 , the accretion flow just above the stellar surface becomes opaque to the stellar 

radiation under the assumption of perfect spherical symmetry. The photosphere is located far from the 
stellar surface, its temperature is reduced to ~ 6000 K even for very massive stars of M* ~ 100 M© (11). 
With the realistic disk accretion, however, the stellar surface at high latitude is not totally embedded in 
such an opaque flow and the UV radiation can freely radiate. For this reason, we evaluate the protostellar 
EUV and FUV photon number luminosities using the fluxes at the stellar surface as 



Seuv = ^Rt I nBi J cs) du, (25) 

Jl3.6cV hv 

S FVV = AnRl r 6 ° V du, (26) 

J11.2eV hV 



where i?* is the stellar radius, B(T c s) is the Planck function, and T e g is defined as 

T _(_ T \ 1/4 



T -ni^f; ■ (27) 

where is the stellar luminosity, L acc = GM*M*/R* is the accretion luminosity, and a is Stcfan- 
Boltzmann constant. 

The stellar UV luminosities given by equations (|25|) - (|27|) assume that all the gas reaching the stellar 
surface releases its gravitational energy at the stellar surface. In reality, some fraction of the gravitational 
energy would be radiated away from the circumstellar disk with a lower effective radiation temperature 
than that of the stellar surface. Because the accretion luminosity L acc is only a few 10 % of the stellar 
luminosity L„, when the protostellar feedback begins to influence the mass accretion (see the main article), 
we do not expect a significant change in the overall accretion process. 

The stellar EUV photons are consumed mostly by photoionizing the recombined hydrogen atoms 
within the HII region. As the recombination rate is proportional to the square of density, this consumption 
is most significant in the vicinity of the protostar. In our calculations, however, the flow structure very 
near the protostar is masked by our assumption of a sink cell. We evaluate the EUV consumption rate 
within the sink cell as follows. First, we consider the spherical "evacuation zone", whose radius i? e v is 
smaller than the gravitational radius for the ionized gas, 



and larger than the size of the sink cell, ~ 10 AU. In our calculations the temperature of ionized gas Xhii 
is ~ 3 x 10 4 K just after formation of the HII region and rises slightly as the stellar mass increases. We 
adopt R cv = 30 AU as a fiducial value. The density distribution within the evacuation zone should be well 
approximated as the free-fall flow whose radial density structure follows p(r) oc r~ 3 / 2 . The consumption 
rate of EUV photons within this zone is analytically written as, 

W,ev = ^ r lnf^, (29) 



87r/i 2 GM* V R 

where a is the total recombination rate, fj, = (1 + 4yHe)w p with the helium abundance j/He and proton 
mass m p , and M ev is the mass inflow rate into this zone. We evaluate M cv using the gas inflow rate along 
the Z-axis, 

32 /2GM7 



47rp z , ov i?^/-^, (30) 



17 



where p z , C v is the gas density at (R,Z) = (0,R CV ). We suppose that the HII region is quenched inside 
the evacuation zone when the stellar EUV luminosity Seuv is less than the consumption rate Seuv.cv 
We do not solve the EUV radiative transfer for such cases. After S'euv exceeds S'euv, ev, we remove the 
limit of Seuv.cv and assume that the EUV luminosity from the sink is Seuv- 

We only marginally resolve the evacuation zone with the current grid resolution. The estimated value 
of Seuv.cv depends on the grid size without resolving i? g .Hii, because free-fall flow is valid only inside 
of Rg.mi- For test calculations with a 2x coarser innermost grid around the protostar, formation of the 
HII region begins for a bit higher stellar mass (by a few 10%). Moreover, the flow structure within the 
evacuation zone could be complex. For instance, the protostellar outflow might be launched from the 
innermost part of the disk (50) with the help of magnetic field generated by dynamo amplification (51, 52). 
This outflow would help the breakout of the HII region by clearing out materials close to the star in the 
polar directions. However, the evolution should depend on the detailed density structure in the outflow- 
launching region, which controls the EUV consumption rate. 

2 Simulation Setup 

As the initial condition of our calculation, we assume the structure of a dense core in the run-away 
collapse from the cosmological simulation by (6). The calculation by (6) followed the entire evolution 
from the cosmological initial condition to the birth of a primordial protostar under the standard ACDM 
cosmology. A small protostar of M* ~ 0.01 Mq forms at 10 20 cm -3 as a result of the run-away collapse 
of a dense primordial-gas core at the cosmological redshift z = 14. Specifically, we take the central 0.3 pc 
cube around the density peak when the maximum density is 10 6 cm -3 as our initial condition. We reduce 
the 3D data to an axisymmetric 2D distribution by averaging over azimuthal angles. Our simulation box 
contains the total gas mass of ~ 300 Mq. 

The numbers of the grid cells are initially Nz x Nr = 42 x 42, including 2 ghost cells in each direction. 
After we start our axisymmetric calculation, the dense core experiences continued gravitational collapse. 
We increase the grid resolution for the central collapsing region by successively adding finer nested grids 
as needed. With this nested-grid technique, we always resolve the minimum Jeans length by tens of 
grid cells. The increase of the grid-level is limited up to 8 owing to computational cost for following 
the subsequent accretion phase until the final stellar mass is fixed. We terminate the calculation of the 
run-away collapse when the minimum Jeans length becomes too short to be resolved with the finest grids. 
At this point we create a sink cell at the origin and calculate the subsequent accretion phase as described 
in Sec. 1. The central density at this moment is ~ 10 12 cm -3 , and the cell size at the finest grid-level is 



This recipe enables us to smoothly connect the 3D cosmological simulation using the particle method 
to a 2D (axial symmetry assumed) local simulation using the nested-grid method. We have confirmed 
that the evolution in our 2D calculation is reasonably consistent with the 3D results after the maximum 
density exceeds 10 6 cm~ 3 . For example, the upper panel of Figure [S3] shows a comparison of the angular- 
momentum profiles for the 3D and 2D calculations using a parameter, 



where M r is the enclosed mass within radius r, and l r is specific angular momentum averaged over the 
spherical shell whose radius is r. The profile for the 2D simulation is within 20 % offset from that in the 
3D case taken from (6). 

Although most of our simulations were done with the settings described above, we also calculated 
early evolution, turning off the stellar feedback, with the higher maximum grid- level of 9. The spatial 
resolution for R < 240 AU is doubled in this case and size of the finest cell is ~ 6 AU. Figure [Si] shows the 



~ 12 AU. 




V<p,r 
Vkep.r 




(31) 
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evolution of the accretion rate until the stellar mass reaches 30 Mq. The accretion rates are initially a 
bit lower with the higher central resolution but converge to our fiducial case as the stellar mass increases. 



3 Structure of the Circumstellar Disk 

Here, we examine the detailed structure of the circumstellar disk observed in our fiducial case. Figure [S4l 
shows the 2D structure of the region within 1000 AU of the protostar, when the stellar mass is ~ 10 M Q 
and 20 Mq. The snapshot for the 10 Mq star shows the structure before the formation of an HII region. 
We see that the accreting gas hits the disk surface and heats up to ~ 4000 K by compression. Radiative 
cooling via H~ free-bound emission operates here as explained in Sec. I1.1.T1 Figure [S5l displays the radial 
structure of the disk at this moment. The hydrogen is primarily in molecular form within the disk due to 
rapid formation via the gas phase three-body reaction. Note that the existence of H2 molecules and their 
contribution for radiative cooling are not included in the semi-analytic models. The total mass of this 
molecular gas is ~ 10 Mq, which is comparable to the stellar mass. Gas on the equator at R < 1000 AU 
has the rotational velocity more than 80% of the local Kepler value. The disk scale height is resolved 
on average except for the innermost part at R < 100 AU. The profile of the Toomre Q-parameter, a 
measure of the gravitational stability of the disk, shows that the minimum value Q m i n is almost unity. 
Thus, the disk is marginally gravitationally stable against the fragmentation. However, high-resolution 
3D simulations demonstrate that the disk fragments due to gravitational instability in very early stage 
of the accretion phase (28). Note that "ring-like" fragmentation can be captured even in axisymmetric 
2D simulations. The absence of fragmentation in our simulation is probably due to the limited spatial 
resolution and the large values of the viscous a-parameter, ag > 0.3 (also see Sec. 11.1.3]) . In our case, 
the disk is also marginally stable in an earlier stage when the stellar mass is ~ 3 Mq (Fig. IS5j) . This is 
almost the same even with the doubled spatial resolution for R < 240 AU. Even disks which experience 
fragmentation have a quasi-static structure with Q m - ln ~ 1, if averaged over many rotation periods. Such 
structure is well mimicked in our calculations. 

The right panel in Figure [S4] shows the 2D disk structure just after the birth of the HII region. 
We see a strong photoevaporating flow in the polar directions within the HII region. The stellar FUV 
luminosity has also increased significantly, but the H2 molecular disk still exists, shadowing the stellar 
photodissociating photons from the dense equatorial regions, evident in the disk's temperature and density 
distributions. The mass of the molecular gas is ~ 10 Mq at this moment. 

4 Comparison to Semi-analytic Models 

Our calculations show that mass accretion toward forming first stars is shut off via dynamical expansion 
of the HII region and photoevaporation of the circumstellar disk. This is qualitatively consistent with 
the picture predicted by the semi-analytic models (8). Their models show that the final stellar mass is 
~ 145 Mq for their fiducial case; this value varies with different sets of input parameters. One of the 
input parameters is /Kep given by equation (|3ip; their adopted fiducial value is 0.5. Figure [S3] shows 
that /kcp — 0.6 — 0.7 for our fiducial case and is lower than 0.5 for the weak-rotation case. For the 
semi-analytic model with /kop > 0.25, however, the final stellar mass does not depend strongly on /kepi 
but rather on the entropy of the accreting gas, measured with a non-dimensional parameter 



where 7 = 1.1 is adopted for a typical value for the primordial gas (2). 

The semi-analytic models use K' = 1 for their fiducial case and show that the final mass roughly scales 
as K' 13 . From Figure [S3l we see that K' ~ 0.7 for our fiducial case and K' ~ 0.8 — 1 for the weak-rotation 




(32) 
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case. The semi-analytic model with K' — 0.7 predicts a final mass ~ 90 M , twice our value of 45 Mq. 
The final mass ~ 85 Mq for our weak-rotation case is somewhat lower than the semi-analytic prediction 
> 108 M Q for K' > 0.8. 

The initial conditions in our examined cases correspond to typical values for gas clouds bearing the 
first stars. Figure [Si] shows that the accretion rates for the cases considered lie in the range 10~ 2 to 
10~ 3 Mq yr _1 for the 10 Mq star, comparable to that expected for 12 mini-halos using the semi-analytic 
models (53). However, our calculations attain systematically lower final masses than predicted by semi- 
analytic models. 

We examine possible reasons for this difference below. First, in our calculations, the mass ratio 
between the star and disk fd = Mdisk/M* differs from the fiducial value in the semi-analytic models, 
fd = 0.3. As discussed in Sec. [3J for our fiducial case the mass of the molecular disk is ~ 10 Mq, 
when the stellar mass is 10 Mq, and would be even higher, if we include the outer atomic part. In 2D 
simulations, however, the value of fd depends on the degree of a- viscosity, an (see Sec. I1.1.3[> . With 
the larger value ao = 1, for example, the mass of the molecular disk for the 10 Mq star is reduced to 
8 Mq. To derive the appropriate value of fd, we need detailed 3D numerical simulations which solve the 
transport of angular momentum in self-gravitating disks. Recent work on the Galactic star formation 
demonstrates that massive disks with fd > 1 can form in some situations, though the value of fd varies 
with different initial conditions of gas cores (54, 55). 

We contend, however, that differences in the stellar feedback process are principally responsible for 
the differences in final stellar mass. The semi-analytic models adopt the following formula for estimating 
the photoevaporation rate of the disk, 

where M„ d — M* + M disk (56). For our fiducial case the photoevaporation rate is about 2 ~ 3 x 
10 -4 Mq yr _1 , after the stellar mass exceeds 30 Mq. By contrast, equation (|3"3")l estimates an evaporation 
rate of several 10 -5 M Q yr" 1 , using Thii = 5 x 10 4 K and the values of S'euv and M*d from the simulation 
(Mdisk is the mass of the molecular disk). We attribute this difference to the basic assumption of an 
infinitely thin disk for deriving equation (f33|) (56). As a result, only the diffuse EUV radiation (see 
Sec. Il.l.l|) is accounted for in equation (p3| . In our calculations the stellar direct EUV radiation field 
impinging on the flared disk is primarily responsible for determining the photoevaporation rate. Even if 
the transfer of diffuse EUV radiation is turned off in the simulation, we do not see any significant change 
in the evolution. 

Moreover, the semi-analytic models assume that, after the HII region forms, mass accretion onto 
the circumstellar disk still continues from regions shaded by the disk, where the stellar UV radiation 
is blocked. In this picture the stellar final mass is determined when the photoevaporation rate exceeds 
the mass accretion rate from the shaded regions. This does not agree with our simulations. It is true 
that the stellar UV radiation is blocked behind the disk as shown in Figure 2 in the main article. In 
our calculations, however, a shock front propagates into the accretion envelope shaded by the disk, as 
the HII region dynamically expands. Figure [Sol shows that there is the outward pressure gradient within 
the HII region due to the photoevaporating flow. As the HII region dynamically expands, the shocked 
accretion envelope behind the disk also obtains the same pressure gradient. As a result, the shocked gas 
is accelerated outward due to this pressure gradient at a velocity of several km/sec. This is shown in 
Figure IS7I We see that the gas at R > a few 10 3 AU is moving outward when the stellar mass exceeds 
20 Mq. The outflow rate of atomic gas via this horizontal motion is 9 x 10~ 4 Mq yr" 1 , when the stellar 
mass is 35 M Q . The mass supply from the envelope to the disk is cut off by this moment H The isolated 

2 In our simulations, we adopt a semi-permeable outer boundary condition, which allows outflow but prohibits mass inflow 
through the boundary. In principle this boundary condition can ultimately limit mass accretion onto the disk. However, 
mass infall toward the disk is reversed and the gas flows outward through the outer boundary long before this limit is 
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disk without replenishment from the envelope continues to lose gas by the photoevaporation. Ultimately, 
mass accretion onto the protostar ceases and the stellar final mass is fixed. 

5 Potential Effects of 3-Dimensionality 

For this study we have adopted an axisymmetric 2D code for calculating the dynamical evolution of 
the accreting flow. Recent 3D simulations for the present-day high-mass star formation are helpful for 
speculating on potential 3D effects. 

A primary radiative feedback effect from Galactic high-mass stars is radiation pressure exerted on 
dust grains coupled with the gas accretion flow. Historically, this effect was first studied under spherical 
symmetry (57-59) and it was found that the accretion flow toward high-mass stars is prevented by the 
strong repulsive force (called as "radiation pressure problem"). 2D numerical simulations show that 
non-spherical mass accretion via disks significantly alleviates this problem (10, 60). 3D simulations 
by (47) demonstrate that 3D effects (e.g., radiative Ray leigh- Taylor instability growing in the accretion 
envelope) also reduce the repulsive force in addition to the effect of the disk accretion. However, it is still 
controversial whether the 3D effects are essential for determining the upper mass limit for present-day 
high- mass star formation. More recent 3D simulations by (61) show that the radiation pressure barrier 
is circumvented with the disk accretion, but find no critical 3D effects. 

Without dust grains the opacity of the primordial gas is much lower than for the present-day Galactic 
interstellar medium. The dynamical evolution discussed in our article is not primarily due to radiation 
pressure, but rather to the high gas pressure of photoionized gas, which causes the dynamical expansion of 
the HII region supplemented by the photoevaporation of the disk. UV feedback effects in the present-day 
high- mass star forming regions have also been studied with 3D numerical simulations (62, 63). However, 
these studies focus on the feedback effects over the large length-scale of cluster-forming clumps, ~ 1 pc. 
By contrast, we consider here the smaller-scale feedback processes in the vicinity of protostars, such as 
photoevaporation of the disks. The photoevaporation rate from the disk is expected to be highest at a 
radial distance from the star R = i? g ,Hii given by equation ([28| (56). However, the regions smaller than 
1000 AU around forming stars are not spatially resolved in the above-cited studies. The photoevaporation 
of the disk is a local process which takes place in the very central part of the HII region. In 3D we can 
expect the photoionized gas to have a clumpier structure than modeled by our 2D simulations. However, 
such small-scale inhomogeneities would be smoothed out within the sound crossing time, which is much 
shorter than the dynamical timescale of the expanding HII region. 

As referred in Sec. 11.1.31 recent 3D simulations demonstrate that the circumstellar disks around 
primordial protostars fragment in the early stage of the accretion phase. The fragmentation of the 
disk would somewhat reduce the stellar final masses, as the accreting materials are shared among the 
fragments (64). Several authors have also suggested that a turbulent velocity field is present after the 
gravitational collapse of primordial gas cores (3). Although the expected turbulent field is not strong for 
the very first stars, we note that, in our case, the turbulent field is further weaken by mapping the 3D data 
to the 2D axisymmetric data (see Sec. [2]). With a random velocity field the angular momentum directions 
of accreting materials would be time-dependent. As a result, the polar direction of the circumstellar disk 
should vary with time as well and the photoevaporating flow would clear out a larger amount of accreting 
materials than otherwise. We speculate that this effect would accelerate the horizontal expansion of the 
HII region and reduce the resulting final stellar mass. Our obtained final mass of the first stars, several 
10 Mq, is thus a conservative upper limit. Nonetheless, this is still much lower than the postulated high 
values exceeding 100 Mq. 
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Fig. SI: Evolution of accretion rates onto the protostar with different a-parameters for angular momentum 
transport. The blue, red, and green curves depict the results with ao — 1, 0.6, and 0.3 in equation (f20|) . 
respectively. The magenta curves display the evolution with the initial angular momentum reduced to 
30% of the fiducial value. For each case the solid and dashed lines represent the evolution with and 
without radiative feedback from the protostar, respectively. The red filled squares represent the no- 
feedback case with a = 0-6 doubling the spatial resolution in the central region of R < 240 AU (also see 
Sec.©. 
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Fig. S2: Accretion rates onto the protostar as a function of the accreted stellar mass, assuming different 
a-parameters for angular momentum transport. The blue, red, and green curves depict the evolution 
with ao = 1, 0.6, and 0.3 in equation (I2TJ1) . respectively. The magenta curves display the evolution with 
the initial angular momentum reduced to 30% of the fiducial value. For each case the solid and dashed 
lines represent the evolution with and without radiative feedback from the protostar, respectively. 
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Fig. S3: Structure of the accreting envelope at the end of the calculation of the run- away collapse stage 
for the fiducial case (red lines) and weak- rotation case (magenta lines). Upper panel: the ratio of the 
local rotational velocity to Kepler velocity /kop- The red open circles represent a snapshot at the birth 
of a embryo protostar in a three-dimensional cosmological simulation (6). Lower panel: Dimensionless 
entropy of the accreting gas K' . The mass- weighted average values K' are plotted against the enclosed 
gas mass. 
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Fig. S4: The structure of gas temperature (left-hand side) and of number density and velocity (right- 
hand side) in the vicinity of 1000 AU around the protostar. The left and right panels show the snapshots 
when the stellar mass is ~ 10 M Q and 20 M©. The elapsed time since the birth of the protostar is also 
shown in each panel. Note that the color scale of the density, the legend of the velocity vectors, and size 
of the plotted area are different from those in Figure 2 in the main article. 
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Fig. S5: Radial structure of the circumstellar disk within 1000 AU around the protostar when the stellar 
mass is 3 Mq (green lines for the lower two panels), 10 Mq (red lines) and 20 Mq (blue lines). Upper 
panel: The fraction of hydrogen molecules along the equator /h 2 (solid lines) , and the ratio of the local 
rotational velocity to Kepler velocity /kep (dashed lines). Middle panel: The scale height of the disk H. 
The black thin line shows the grid resolution in our standard cases, and the discontinuity at R ~ 470 AU 
indicates a grid-level boundary there. The profile when the stellar mass is 3 Mq with the doubled spatial 
resolution for R < 240 AU is also presented (green filled squares). The black thin dashed line shows the 
grid resolution in this case. Lower panel: Same as the middle panel but for Toomre Q-parameter at each 
radius. 
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Fig. S6: Same as Fig. 2 in the main article but for gas temperature multiplied by number density in 
the left panel, which is nearly proportional to gas pressure. The left and right panels show snapshots at 
times, when the stellar mass is M* = 32 M and 37.5 M Q respectively. 
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Fig. S7: Distribution of i?-component of the radial velocity Vr along the equator (Z = 0). The snapshots 
when the stellar mass is 10 M (black), 20 M Q (red), 32 M Q (blue), 37.5 M Q (magenta), and 42 M e 
(green) are plotted. Positive vr means that gas is moving outward. 
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